﻿Imports System.Collections.Generic
Imports System.Collections.ObjectModel
Imports System.Linq
Imports System.Text
Imports DotSpatial.Controls
Imports GeoAPI.Geometries

'using DotSpatial.Data
Imports NetTopologySuite.Geometries
Imports NetTopologySuite.Noding.Snapround
Imports NetTopologySuite.Operation.Linemerge
Imports NetTopologySuite.Operation.Polygonize

' ********************************************************************************************************
' Product Name: Conturer
' Description:  An extension for create contour from raste
' ********************************************************************************************************
' The contents of this file are subject to the MIT License (MIT)
' you may not use this file except in compliance with the License. You may obtain a copy of the License at
' http://dotspatial.codeplex.com/license
'
' Software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTY OF
' ANY KIND, either expressed or implied. See the License for the specific language governing rights and
' limitations under the License.
'
' The Initial Developer of this Original Code is Tomaso Tonelli in memory of his mother in 2012.
' for questions/ info tomaso.tonelli@gmail.com
' ********************************************************************************************************

''' <summary>
''' Converted from C# source code DotSpatial.Plugins.Contourer version 70488
''' </summary>
''' <remarks></remarks>
Public Class clsContour
    Public Enum ContourType
        Line
        Polygon
    End Enum

    Public Shared Function Execute(rst As DotSpatial.Data.Raster, type As ContourType, Optional FieldName As String = "Value", Optional levels As Double() = Nothing) As DotSpatial.Data.FeatureSet
        Dim lev As Double() = levels

        Dim iRst As DotSpatial.Data.Raster = RasterCheck(rst, lev)



        Dim field As String

        If FieldName Is Nothing Then
            field = "Value"
        Else
            field = FieldName
        End If

        Dim x As Double() = New Double(rst.NumColumns - 1) {}
        Dim y As Double() = New Double(rst.NumRows - 1) {}

        For i As Integer = 0 To rst.NumColumns - 1
            x(i) = rst.Extent.MinX + rst.CellWidth * i + rst.CellWidth / 2
        Next

        For i As Integer = 0 To rst.NumRows - 1
            y(i) = rst.Extent.MaxY - rst.CellHeight * i - rst.CellHeight / 2
        Next

        Dim fs As DotSpatial.Data.FeatureSet = Nothing

        Select Case type
            Case ContourType.Line
                If True Then
                    fs = New DotSpatial.Data.FeatureSet(DotSpatial.Topology.FeatureType.Line)
                    fs.DataTable.Columns.Add(field, GetType(Double))

                    For z As Integer = 0 To levels.Count() - 1
                        Dim cont As IList(Of IGeometry) = GetContours(iRst, x, y, lev(z))

                        For Each g As Geometry In cont
                            Dim f As DotSpatial.Data.Feature = DirectCast(fs.AddFeature(ToDotSpatialLineString(DirectCast(g, ILineString))), DotSpatial.Data.Feature)
                            f.DataRow(field) = lev(z)
                        Next
                    Next
                End If
                Exit Select

            Case ContourType.Polygon
                If True Then
                    fs = New DotSpatial.Data.FeatureSet(DotSpatial.Topology.FeatureType.Polygon)

                    fs.DataTable.Columns.Add("Lev", GetType(Integer))
                    fs.DataTable.Columns.Add("Label", GetType(String))

                    Dim Contours As New Collection(Of IGeometry)()

                    For z As Integer = 0 To levels.Count() - 1
                        Dim cont As IList(Of IGeometry) = GetContours(iRst, x, y, lev(z))

                        For Each g As Geometry In cont
                            Contours.Add(New LineString(g.Coordinates))
                        Next
                    Next

                    Dim Boundary As Coordinate() = New Coordinate(4) {}

                    Boundary(0) = New Coordinate(x(0), y(0))
                    Boundary(1) = New Coordinate(x(0), y(rst.NumRows - 1))
                    Boundary(2) = New Coordinate(x(rst.NumColumns - 1), y(rst.NumRows - 1))
                    Boundary(3) = New Coordinate(x(rst.NumColumns - 1), y(0))
                    Boundary(4) = New Coordinate(x(0), y(0))

                    Contours.Add(New LineString(Boundary))

                    Dim NodedContours As New Collection(Of IGeometry)()
                    Dim geomNoder As New GeometryNoder(New PrecisionModel(1000.0))

                    For Each c As LineString In geomNoder.Node(Contours)
                        NodedContours.Add(c)
                    Next

                    Dim polygonizer As New Polygonizer()
                    polygonizer.Add(NodedContours)

                    For Each p As Polygon In polygonizer.GetPolygons()
                        Dim pnt As Point = DirectCast(p.InteriorPoint, Point)

                        Dim c As Integer = CInt((pnt.X - iRst.Extent.MinX) / iRst.CellWidth)
                        Dim r As Integer = CInt((iRst.Extent.MaxY - pnt.Y) / iRst.CellHeight)

                        Dim z As Double = DirectCast(iRst, DotSpatial.Data.Raster).Value(r, c)

                        Dim Cls As Integer = GetLevel(z, lev)
                        Dim label As String = "Undefined"

                        If Cls = -1 Then
                            label = "< " + lev(0).ToString()
                        End If
                        If Cls = lev.Count() Then
                            label = "> " + lev(lev.Count() - 1).ToString()
                        End If
                        If Cls >= 0 And Cls < lev.Count() Then
                            label = lev(Cls).ToString() + " - " + lev(Cls + 1).ToString()
                        End If

                        Dim dsp As DotSpatial.Topology.Polygon = ToDotSpatialPolygon(p)

                        Dim f As DotSpatial.Data.Feature = DirectCast(fs.AddFeature(dsp), DotSpatial.Data.Feature)
                        f.DataRow("Lev") = Cls
                        f.DataRow("Label") = label
                    Next
                End If
                Exit Select
        End Select

        Return fs
    End Function

    Public Shared Function GetLevel(z As Double, lev As Double()) As Integer
        If z < lev(0) Then
            Return -1
        End If
        If z > lev(lev.Count() - 1) Then
            Return lev.Count()
        End If

        For i As Integer = 0 To lev.Count() - 2
            If z >= lev(i) And z <= lev(i + 1) Then
                Return i
            End If
        Next

        Return -9999
    End Function

    Public Shared Sub CreateMinMaxEvery(r As DotSpatial.Data.Raster, type As ContourType, ByRef MinContour As Double, ByRef MaxContour As Double, ByRef every As Double)
        Dim min As Double = r.Minimum
        Dim max As Double = r.Maximum

        If min = max Then
            min = Math.Floor(min)
            max = Math.Ceiling(max)
            If min = max Then
                max += 1
            End If
        End If

        Dim dz As Double = max - min

        Dim Order As Double = Math.Pow(10, Math.Floor(Math.Log10(dz)))

        If Order = dz Then
            Order /= 10
        End If

        If dz / Order < 2 Then
            Order /= 10
        End If

        MinContour = Math.Floor(min / Order) * Order
        MaxContour = Math.Ceiling(max / Order) * Order

        If MaxContour < max Then
            MaxContour += Order
        End If

        every = Order

        If type = ContourType.Line Then
            MinContour += every
            MaxContour -= every

            If MaxContour <= MinContour Then
                MaxContour = MinContour + every
            End If
        End If
    End Sub

    Public Shared Function CreateLevels(MinContour As Double, MaxContour As Double, every As Double) As Double()
        Dim c As Integer = CInt((MaxContour - MinContour) / every) + 1

        Dim levels As Double() = New Double(c - 1) {}

        For i As Integer = 0 To c - 1
            levels(i) = MinContour + every * i
        Next

        Return levels
    End Function

    Public Shared Function toDotSpatial(l As Geometry) As DotSpatial.Topology.Geometry
        If l.GeometryType = "LineString" Then
            Return ToDotSpatialLineString(DirectCast(l, ILineString))
        End If

        If l.GeometryType = "Polygon" Then
            Dim p As Polygon = TryCast(l, Polygon)

            Dim dsp As DotSpatial.Topology.Polygon = ToDotSpatialPolygon(DirectCast(l, GeoAPI.Geometries.IPolygon))
            Return DirectCast(dsp, DotSpatial.Topology.Geometry)
        End If

        Return Nothing
    End Function

    Friend Shared Function ToDotSpatialPoint(point As GeoAPI.Geometries.IPoint) As DotSpatial.Topology.Point
        Return New DotSpatial.Topology.Point(point.X, point.Y)
    End Function

    Friend Shared Function ToDotSpatialPoint(point As GeoAPI.Geometries.Coordinate) As DotSpatial.Topology.Point
        Return New DotSpatial.Topology.Point(point.X, point.Y)
    End Function

    Friend Shared Function ToDotSpatialCoordinate(coordinate As GeoAPI.Geometries.Coordinate) As DotSpatial.Topology.Coordinate
        Return New DotSpatial.Topology.Coordinate(coordinate.X, coordinate.Y)
    End Function

    Friend Shared Function ToDotSpatialLineString(l As GeoAPI.Geometries.ILineString) As DotSpatial.Topology.LineString
        Dim c As DotSpatial.Topology.Coordinate() = New DotSpatial.Topology.Coordinate(l.Coordinates.Count() - 1) {}

        For i As Integer = 0 To l.Coordinates.Count() - 1
            c(i) = New DotSpatial.Topology.Coordinate(l.Coordinates(i).X, l.Coordinates(i).Y)
        Next
        Return New DotSpatial.Topology.LineString(c)
    End Function

    Friend Shared Function ToDotSpatialLinearRing(geom As GeoAPI.Geometries.ILinearRing) As DotSpatial.Topology.LinearRing
        Dim vertices As New Collection(Of DotSpatial.Topology.Point)()

        For Each coordinate As Coordinate In geom.Coordinates
            Dim p As DotSpatial.Topology.Point = ToDotSpatialPoint(coordinate)

            vertices.Add(p)
        Next
        Return New DotSpatial.Topology.LinearRing(vertices)
    End Function

    Friend Shared Function ToDotSpatialPolygon(geom As GeoAPI.Geometries.IPolygon) As DotSpatial.Topology.Polygon
        Dim exteriorRing As DotSpatial.Topology.LinearRing = ToDotSpatialLinearRing(DirectCast(geom.ExteriorRing, GeoAPI.Geometries.ILinearRing))

        Dim interiorRings As DotSpatial.Topology.LinearRing() = New DotSpatial.Topology.LinearRing(geom.InteriorRings.Count() - 1) {}

        'foreach (GeoAPI.Geometries.ILineString interiorRing in geom.InteriorRings)
        For i As Integer = 0 To geom.InteriorRings.Count() - 1
            Dim hole As DotSpatial.Topology.LinearRing = ToDotSpatialLinearRing(DirectCast(geom.InteriorRings(i), GeoAPI.Geometries.ILinearRing))
            interiorRings(i) = hole
        Next

        Return New DotSpatial.Topology.Polygon(exteriorRing, interiorRings)
    End Function

    Public Shared Function RasterCheck(raster As DotSpatial.Data.Raster, levels As Double()) As DotSpatial.Data.Raster
        Dim eps As Double = (raster.Maximum - raster.Minimum) / 1000
        Dim n As Integer = levels.Count()

        Dim rst As DotSpatial.Data.Raster = raster

        For i As Integer = 0 To rst.NumRows - 1
            For j As Integer = 0 To rst.NumColumns - 1
                If rst.Value(i, j) <> rst.NoDataValue Then
                    For l As Integer = 0 To n - 1
                        If rst.Value(i, j) = levels(l) Then
                            rst.Value(i, j) += eps
                        End If
                    Next
                End If
            Next
        Next

        Return rst
    End Function

    Public Shared Function GetContours(ByRef rst As DotSpatial.Data.Raster, x As Double(), y As Double(), zlev As Double) As IList(Of IGeometry)
        Dim lsList As New List(Of LineString)()

        Dim ipari As Boolean, jpari As Boolean
        Dim xx As Double() = New Double(2) {}
        Dim yy As Double() = New Double(2) {}
        Dim zz As Double() = New Double(2) {}

        For j As Integer = 0 To rst.NumColumns - 2
            If CInt(CDbl(j) / 2.0) * 2 = j Then
                jpari = True
            Else
                jpari = False
            End If

            For i As Integer = 0 To rst.NumRows - 2
                If CInt(CDbl(i) / 2.0) * 2 = i Then
                    ipari = True
                Else
                    ipari = False
                End If

                If Not jpari AndAlso Not ipari OrElse jpari AndAlso ipari Then
                    xx(0) = x(j)
                    yy(0) = y(i)
                    zz(0) = rst.Value(i, j)

                    xx(1) = x(j)
                    yy(1) = y(i + 1)
                    zz(1) = rst.Value(i + 1, j)

                    xx(2) = x(j + 1)
                    yy(2) = y(i)
                    zz(2) = rst.Value(i, j + 1)

                    If True Then
                        Dim c As Coordinate() = Intersect(xx, yy, zz, zlev)
                        If c IsNot Nothing Then
                            lsList.Add(New LineString(c))
                        End If
                    End If

                    xx(0) = x(j + 1)
                    yy(0) = y(i)
                    zz(0) = rst.Value(i, j + 1)

                    xx(1) = x(j)
                    yy(1) = y(i + 1)
                    zz(1) = rst.Value(i + 1, j)

                    xx(2) = x(j + 1)
                    yy(2) = y(i + 1)
                    zz(2) = rst.Value(i + 1, j + 1)

                    If True Then
                        Dim c As Coordinate() = Intersect(xx, yy, zz, zlev)
                        If c IsNot Nothing Then
                            lsList.Add(New LineString(c))
                        End If
                    End If
                End If

                If jpari AndAlso Not ipari OrElse Not jpari AndAlso ipari Then
                    xx(0) = x(j)
                    yy(0) = y(i)
                    zz(0) = rst.Value(i, j)

                    xx(1) = x(j)
                    yy(1) = y(i + 1)
                    zz(1) = rst.Value(i + 1, j)

                    xx(2) = x(j + 1)
                    yy(2) = y(i + 1)
                    zz(2) = rst.Value(i + 1, j + 1)

                    If True Then
                        Dim c As Coordinate() = Intersect(xx, yy, zz, zlev)
                        If c IsNot Nothing Then
                            lsList.Add(New LineString(c))
                        End If
                    End If

                    xx(0) = x(j)
                    yy(0) = y(i)
                    zz(0) = rst.Value(i, j)

                    xx(1) = x(j + 1)
                    yy(1) = y(i + 1)
                    zz(1) = rst.Value(i + 1, j + 1)

                    xx(2) = x(j + 1)
                    yy(2) = y(i)
                    zz(2) = rst.Value(i, j + 1)

                    If True Then
                        Dim c As Coordinate() = Intersect(xx, yy, zz, zlev)
                        If c IsNot Nothing Then
                            lsList.Add(New LineString(c))
                        End If
                    End If
                End If
            Next
        Next

        'Comapct(lsList)
        'Return lsList

        'look for duplicates
        Dim numdup As Integer = 0
        For i = 0 To lsList.Count - 1
            For j = i + 1 To lsList.Count - 1
                If lsList(i).EqualsExact(lsList(j)) Then numdup += 1
            Next
        Next

        Dim lm As New LineMerger()

        lm.Add(lsList)

        Dim merged As IList(Of IGeometry) = DirectCast(lm.GetMergedLineStrings(), IList(Of IGeometry))

        Return merged
    End Function

    Private Shared Sub Comapct(ByRef lsList As List(Of LineString))
        Dim n As Integer = lsList.Count()
        Dim found As Boolean = True

        While found
            found = False

            For i As Integer = 0 To n - 3
                For j As Integer = i + 1 To n - 2
                    Dim l As LineString = UnionLineStrings(lsList(i), lsList(j))

                    If l IsNot Nothing Then
                        lsList(i) = l
                        lsList.RemoveAt(j)
                        n -= 1
                        found = True
                        If j > n - 3 Then Exit For
                    End If
                Next

            Next
        End While

    End Sub

    Private Shared Function UnionLineStrings(ls1 As LineString, ls2 As LineString) As LineString
        Dim l As LineString = Nothing

        If ls1.Coordinates.Count() < 0 Or ls2.Coordinates.Count() < 0 Then
            Return Nothing
        End If

        Dim c1 As Coordinate() = Nothing
        Dim c2 As Coordinate() = Nothing

        If ls1.StartPoint.Distance(ls2.StartPoint) < 0.001 Then
            c1 = Reverse(ls1.Coordinates)
            c2 = ls2.Coordinates
        End If

        If ls1.EndPoint.Distance(ls2.EndPoint) < 0.001 Then
            c1 = ls1.Coordinates
            'c.Add((CoordinateList)ls2.Coordinates.Reverse(), false);
            c2 = Reverse(ls2.Coordinates)
        End If

        If ls1.StartPoint.Distance(ls2.EndPoint) < 0.001 Then
            c1 = ls2.Coordinates
            'c.Add((CoordinateList)ls1.Coordinates, false);
            c2 = ls1.Coordinates
        End If

        If ls1.EndPoint.Distance(ls2.StartPoint) < 0.001 Then
            c1 = ls1.Coordinates
            'c.Add((CoordinateList)ls2.Coordinates, false);
            c2 = ls2.Coordinates
        End If

        If c1 IsNot Nothing AndAlso c2 IsNot Nothing Then
            Dim xx As Coordinate() = c1
            xx.Concat(c2)
            l = New LineString(xx)
        End If
        Return l
    End Function

    Private Shared Function Reverse(l As Coordinate()) As Coordinate()
        Dim x As Coordinate() = l

        Array.Reverse(x)

        Return x

    End Function

    Private Shared Function Intersect(xx As Double(), yy As Double(), zz As Double(), zlevel As Double) As Coordinate()
        Dim coordinates As New List(Of Coordinate)()

        Dim zmin As Double = Math.Min(zz(0), Math.Min(zz(1), zz(2)))
        Dim zmax As Double = Math.Max(zz(0), Math.Max(zz(1), zz(2)))

        If zlevel < zmin OrElse zlevel > zmax Then
            Return Nothing
        End If

        Dim i As Integer, i1 As Integer

        For i = 0 To 2
            i1 = i + 1
            If i = 2 Then
                i1 = 0
            End If

            zmin = Math.Min(zz(i), zz(i1))
            zmax = Math.Max(zz(i), zz(i1))

            If zlevel > zmin AndAlso zlevel <= zmax Then
                Dim dz As Double = zz(i1) - zz(i)
                If True Then
                    Dim frac As Double = (zlevel - zz(i)) / dz

                    If frac > 0 Then
                        Dim dx As Double = (xx(i1) - xx(i))
                        Dim dy As Double = (yy(i1) - yy(i))

                        Dim c As New Coordinate()

                        c.X = xx(i) + frac * dx
                        c.Y = yy(i) + frac * dy
                        c.Z = zlevel

                        coordinates.Add(c)
                    End If
                End If
            End If
        Next

        If coordinates.Count >= 3 Then
            Throw New Exception("three intersection!")
        End If

        If coordinates.Count = 1 Then
            Throw New Exception("one intersection!")
        End If

        Return coordinates.ToArray()
    End Function
End Class
